/*
 *  mooneyrivlin.cpp
 *  EPPI-0.0
 *
 *  Created by Sergio Conde Martin on 14/02/2013.
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */

#include "mooneyrivlin.h"

MooneyRivlin::MooneyRivlin()
:	HyperelasticInvariants(),
 	c1(0.0),
 	c2(0.0),
 	b(0.0),
	typeVol("")
{}

// ***

MooneyRivlin::MooneyRivlin(ifstream &in, const string& na, const string & t)
:	HyperelasticInvariants(na,t),
 	c1(0.0),
 	c2(0.0),
 	b(0.0),
	typeVol("")
{
	string keyword;
	while (in >> keyword && keyword.compare("ENDMATERIAL") )
	{
		if(!keyword.compare("VOL"))
		{
			in >> typeVol; 
			if(!typeVol.compare("mixed")) in >> b;
		}
		else if(!keyword.compare("C1"))	in >> c1;
		else if(!keyword.compare("C2"))	in >> c2;
        else if(!keyword.compare("K"))  in >> K;
		else if(!keyword.compare("DENSITY"))	in >> density;
	}

	mu = 2.0 * (c1 + c2);
    E = 9.0 * K * mu / (3.0 * K + mu);
    nu = (3.0 * K -2.0 * mu ) / (2.0 * (3.0 * K + mu));
    lambda = nu * E / ((1.0 + nu) * (1.0 - 2.0 * nu));
	
	hypertype ="coupled";
	
	if(!typeVol.compare("cuad") || !typeVol.compare("log") || !typeVol.compare("mixed"))
	{
		hypertype = "decoupled";
	}
	if (!hypertype.compare("decoupled"))
	{
		Svol.setZero(6);
		Siso.setZero(6);
		Ctvol.setZero(6,6);
		Ctiso.setZero(6,6);
	}
}

// *** 

void MooneyRivlin::strainEnergy()
{
      // Compressible or nearly incompressible mooneyrivlin material (volumetric part as neohooke in Bonet's book -pag. 16-).
	if (J > 0.0) 
	{
		if(!hypertype.compare("coupled"))
		{
			W = c1 * (I - 3.0) + c2 * (II - 3.0) - mu * log(J) + lambda / 2.0 * log(J) * log(J); 
		}
		else
		{
			Wiso = c1 * (Im - 3.0) + c2 * (IIm - 3.0);
			
			if(!typeVol.compare("cuad"))	Wvol = K / 2.0 * (J - 1.0) * (J - 1.0);
			if(!typeVol.compare("log"))     Wvol = K / 2.0 * log(J) * log(J);
			if(!typeVol.compare("mixed"))	Wvol = K / (b * b) * (1.0 / pow(J,b) - 1.0 + b *log(J));
			
			W = Wvol + Wiso;
		}
	}
	else
	{
		cout << "*ERROR: Negative volume, material " << name << " has collapsed \n";
		Wvol = -0.5;
		Wiso = -0.5;
		W = -1.0;
	}
}

// ***

void MooneyRivlin::DstrainEnergy()
{
	if (J > 0.0)
	{
		if(!hypertype.compare("coupled"))
		{
			W1 = c1;
			W2 = c2;
			W3 = ( lambda * log(III) - 2.0 * mu) / (4.0 * III);
		}
		else
		{
			Wiso1 = c1;
			Wiso2 = c2;
			
			if(!typeVol.compare("cuad"))	dWvol = K * (J - 1.0);
			if(!typeVol.compare("log"))     dWvol = K * log(J) / J;
			if(!typeVol.compare("mixed"))	dWvol = K / (b * J) * (1.0 - 1.0 / pow(J,b));
		}
	}
	else
	{
		W1 = -1.0;
		W2 = -1.0;
		W3 = -1.0;
		
		dWvol = -1.0;
		Wiso1 = -1.0;
		Wiso2 = -1.0;
	}
}

// ***

void MooneyRivlin::DDstrainEnergy()
{
	if (J > 0.0)
	{
		if(!hypertype.compare("coupled"))
		{
			W11 = 0.0;
			W12 = 0.0;
			W13 = 0.0;
			W22 = 0.0;
			W23 = 0.0;
			W33 = ( 2.0 * mu + lambda * (1.0 - log(III)) ) / (4.0 * III * III);
		}
		else
		{
			Wiso11 = 0.0;
			Wiso12 = 0.0;
			Wiso22 = 0.0;
			
			if(!typeVol.compare("cuad"))	ddWvol = K; 
			if(!typeVol.compare("log"))     ddWvol = K * (1.0 - log(J)) / (J * J);
			if(!typeVol.compare("mixed"))	ddWvol = K / (b * J * J) * ( (b + 1.0) / pow(J,b) - 1.0);
		}
	}
	else
	{
		W11 = -1.0;
		W12 = -1.0;
		W13 = -1.0;
		W22 = -1.0;
		W23 = -1.0;
		W33 = -1.0;
		
		ddWvol = - 1.0; 
		Wiso11 = - 1.0;
		Wiso12 = - 1.0;
		Wiso22 = - 1.0;
	}
}

// ***

void MooneyRivlin::DDDstrainEnergy()
{}

// ***

//double MooneyRivlin::getE()
//{	return E;}
//
//// ***
//
//double MooneyRivlin::getNu()
//{	return nu;}
//
//// ***
//
//double MooneyRivlin::getLambda()
//{	return lambda;}
//
//// ***
//
//double MooneyRivlin::getMu()
//{	return mu;}
//
//// ***
//
//double MooneyRivlin::getDensity()
//{	return density;}
//
//// ***
//
//double MooneyRivlin::getBulkmodulus()
//{	return K;}
